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1. The Lagrangian Model 

Under this grant we have developed and tested a new Lagrangian model of sea ice. A Lagrangian 
model keeps track of material parcels as they drift in the model domain. Besides providing a natural 
framework for the assimilation of Lagrangian data, it has other advantages: 1) a model that follows 
material elements is well suited for a medium such as sea ice in which an element retains its identity for 
a long period of time; 2) model cells can be added or dropped as needed, allowing the spatial resolution 
to be increased in areas of high variability or dense observations; 3) ice from particular regions, such as 
the marginal seas, can be marked and traced for a long time; and 4) slip lines in the ice motion are 
accommodated more naturally because there is no internal grid. Our work makes use of these strengths 
of the Lagrangian formulation. 

The Lagrangian model is a powerful new tool that could be used to address important science ques- 
tions regarding the movement of sea ice. Ice is known to grow predominantly in selected shelf regions 
where there is significant offshore flow and ice divergence in the winter, and to melt either in shelf 
regions or in the Greenland Sea. How does the ice move from source regions to melt regions and to 
what extent does the ice from a source region mix with multiyear ice from the central pack? How do 
these trajectories change from year to year? How does the ice thickness of a parcel change over a mul- 
tiple-year trajectory? 

Our model is based on a set of Lagrangian cells. The cells are considered as individual regions of 
ice in the dynamic model, and each cell has an associated position, velocity, ice concentration, and 
mean ice thickness. The cells are initially given positions in a square grid and the areas are assigned to 
fill the space at all times. The initial cells may be more densely assigned in regions where higher reso- 
lution is required. Cells may be removed if they lose all their ice, or merged with a neighbor if they 
move too close together. A cell may also be divided into two cells if it becomes too large or if higher 
resolution is required. Cells are created when ice forms in formerly ice-free regions and may also be 
created to correspond to the initial position of observed trajectories. The position and velocity of each 
cell is integrated forward in time with an adaptive time stepping procedure which accounts for the time- 
varying forces acting on the cell. The forces include wind stress, water stress, Coriolis force, and inter- 
nal ice stress. In addition a coastal force is imposed to keep the cells away from land and a corrective 
force may be imposed to force a cell to follow an observed trajectory. 

The spatial resolution of the model is established by the local number density of the cells and ranges 
from 10 km where there are many observations to 100 km where there are none. The basic temporal 
output is 1/2 day, but an adaptive time-stepping procedure breaks the basic step down into 100 to 2000 
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smaller steps in order to obtain the required accuracy. This allows the determination of cell positions to 
within a few minutes of any observed trajectory positions. 

A number of improvements have already been incorporated into the model: The length scale for 
determining the strain rate is now variable and depends on the number density of the cells. Where the 
density is high the length scale is small. The forcing ocean currents are now time varying, not seasonal 
climatological values. The daily averaged currents are taken from the output of the Polar Science Cen- 
ter Eulerian coupled ice/ocean model run with assimilation of buoy and SSMI-based ice displacement 
data (Zhang et al., 2003). This Eulerian model run with velocity data assimilation gives very accurate 
surface currents. Furthermore, the data assimilation procedures outlined below have now been success- 
fully implemented. 

2. Data assimilation procedures 

An important new tool for observing and understanding Arctic sea ice movement and deformation 
has been implemented in the last few years. The RADARSAT Geophysical Processor System (RGPS) 
uses synthetic aperture radar (SAR) images of the Arctic Ocean from the Canadian RADARSAT satel- 
lite as input to an ice-tracking algorithm that identifies and follows more than 40,000 points on the sea 
ice by means of a cross-correlation technique (Kwok 1998). RADARSAT swaths are 460 km wide, and 
complete coverage of the western Arctic Ocean takes about three days to acquire. The ice tracking 
begins each November after freeze-up, with points initially placed on a regular 10-km grid, and contin- 
ues until the following spring, when the onset of melt erases the distinct signatures of the surface fea- 
tures. Each point is tracked in a Lagrangian fashion for an entire season. 

However the sampling of the RGPS is very idiosyncratic. It is irregular in both time and space and 
the temporal sampling, typically three days, is longer than is commonly used in model simulations. The 
irregular sampling makes it very difficult to analyze the deformation fields to determine the spatial 
structure of the deformation patterns. One solution is to average the data spatially and temporally to 
regularize the data set. A monthly averaged gridded data set produced under this grant is posted on the 
Arctic Ocean Model Intercomparison (AOMIP) web site and at psc.apl.washington.edu/lindsay/#RGPS. 
But in this case the fine spatial and temporal resolution is lost. One of our goals is to regularize the 
data set without losing the spatial and temporal resolution. 

In order to obtain the best estimates of the motion of the ice pack with regular time and space sam- 
pling from the RGPS trajectory data, we use novel data assimilation procedures to constrain the state of 
the motion. Specific cells in the Lagrangian model are associated with individual RGPS trajectories and 
will, over the course of the year, be constrained to move as the observations indicate. The assimilation 
is accomplished in a three-part process. 

The first part is to run the model without assimilation for an interval of time, say 20 days. Any 
observed ice trajectory that is to be used, either for assimilation or for validation, is assigned a model 
cell when the trajectory is first available. We often run the assimilation procedures for the even num- 
bered cells and use the odd numbered cells for evaluation. These simulated trajectories reflect the 
model-estimated wind stress, water stress, Coriolis force, and internal stress terms of the momentum 
balance equation. 

The second step constructs a target path for each assimilation cell by adjusting the first-guess trajec- 
tory so that it passes through the observed positions. The adjustment is made by subtracting the linearly 
interpolated position error in the trajectory from the first guess. This preserves some of the shape of the 
first-guess trajectory that is indicative of the changing forces acting on the cell while still forcing the 
target trajectory through the observations. 

The final step is to integrate the interval again and include a corrective force. This force represents 
the net error in various terms of the momentum balance equation, such as the wind and ocean drag or 
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the internal stresses. If ^target > s the target trajectory position (x component) and x mod is the current 
model position, then the corrective force is 
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where k is a positive constant. The Gaussian exponential reduces the force when no observation is 
close in time: t is the time to the nearest observation and T is a constant, currently set at 2 days. Thus 
the force is strongest near an observation time or when the position error is large. The influence of the 
force can be made arbitrarily large by increasing k. Consequently a very good match, less than 1 km, is 
accomplished between the model and the observed trajectories. This corrective force is then interpo- 
lated smoothly to all the Lagrangian cells lacking an accompanying displacement measurement. Evalu- 
ation of the spatial and temporal variation of this corrective force is a key element of our model 
validation efforts. 

The corrective force is perhaps the most interesting aspect of this data assimilation technique. Fig- 
ure 1 shows a short time series of a trajectory (in one dimension) of a single cell: as observed, the first- 
guess, the target, and as finally modeled. Four data assimilation periods are shown. The first-guess tra- 
jectory (red), with no data assimilation, wanders away from the observed trajectory, although there is 
some correspondence. The second guess trajectory (blue) is almost overlain by the final model result 
(green). The final model result passes through the observed positions of the trajectory. The lower panel 
shows the five principal forces acting on the cell. The corrective force (black) is largest when the sec- 
ond-guess and final result differ. At about day 326 the corrective force is large because the wind stress 
(cyan) is not strong enough to push the cell south as rapidly as observed. In contrast, on day 339 the 
corrective force is in the opposite direction from the large positive wind stress and it is needed to coun- 
teract the wind stress because the cell moves very little during this period. Uncertainties in the other 
forces, including the internal stress, are also incorporated in the corrective force so that errors in the 
wind stress may not be the only contributor to the corrective force. In our preliminary analysis of the 
corrective force, we have been struck by the fact that the corrective force is largest, and often in the 
opposite sense, when the wind stress is high. 

The corrective force is a measure of the error in the vector sum of the forces mentioned above. Fig- 
ure 2 shows an example of the forces computed for cells within a 500-km region. We see that there is 
spatial coherence in the corrective force field but it is far from constant and exhibits a strong discontinu- 
ity. We believe that the discontinuity arises from a shear crack in the ice, a linear kinematic feature 
(Kwok, 2001), that is characterized by a discontinuity in the velocity. This discontinuity in the velocity 
required a discontinuity in the corrective forces that is likely associated with incorrect internal stress 
gradient estimates in the model. 

A definitive partition of the error among the forces is not possible. However it will be instructive to 
examine the mean field of the corrective force and how it relates to the mean fields of the other forces 
in the model. Figure 3 shows a preliminary estimate of the mean ice velocity and the mean of three of 
the forces. We see that the mean corrective force has a distinctive spatial pattern that is not a simple 
reflection of the wind stress and may be indicating significant errors in the internal stress gradient, since 
it is strongest near the Alaskan coast and where the ice is converging against the Canadian islands. The 
ice strength in this simulation is perhaps too small. 

Validation of the model ice velocity is accomplished by comparing interpolated velocities with daily 
averaged buoy velocities from the LABP. The excellent temporal resolution of the buoy data provide 
good validation for our nominal 1-day sampling interval. For example, from the SHEBA year we have 
produced a 50-km data set by assimilating a random subset of the RGPS trajectories. The RMS error of 
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the trajectories within 10 km of a buoy is 1.0 km/day and the correlation is R = 0.95. Table 1 shows a 
comparison of the error characteristics of different methods of measuring or estimating the ice velocity. 
Note that of these methods the Lagrangian model with assimilation of RGPS trajectories is the best. 


Table 1: Error characteristics of various data sets compared to 24-hour buoy velocities. 



Correlation 

(R) 

RMS error 
(km/day) 

Source 

Buoy measurement error 


0.4 

IABP 

RGPS 3-day displacements 

0.92 

2.3 

Lindsay (2002) 

Eulerian model, no DA 

0.66 

6.0 

Zhang et al., (2003) 

Eulerian Model, w/ DA 
of buoy and SSMI velocities 

0.90 

3.5 

Zhang et al., (2003) 

Lagrangian model, no DA 

0.80 

5.3 

Lindsay and Stem (2004) 

Lagrangian Model, w/ DA of 2500 
RGPS trajectories, SHEBA winter 

0.95 

1.0 

current work 


DA = data assimilation 


3. Accomplishments made under this grant 

• The primary product of the research so far is a paper published in the Journal of Physical Oceanog- 
raphy which is attached as an addendum to this report. This paper describes the new Lagrangian 
model of sea ice and shows how it is able to simulate the ice motion about as well as a state-of-the- 
art Eulerian model (Table 1, lines 3 and 5). A simulation of the ice motion over the five-year period 
1993-1997 is shown. The paper describes the model but does not describe our new assimilation 
techniques. 

• Gave a talk “A Lagrangian Dynamic Model of Sea Ice for Data Assimilation” at the Workshop on 
Sea Ice Data Assimilation Naval Academy, Annapolis MD July 23-24, 2002 

• Helped edit “Report to NASA on a Workshop on Sea Ice Data Assimilation” held at the Naval 
Academy, Annapolis MD, July 2002. 

• Attended the “Short Course on Data Assimilation for Sea-Ice Modelers”, Woods Hole Oceano- 
graphic Institute, Woods Hole, MA, May 10-11, 2003. 

• Presented a poster “A new Lagrangian model of sea ice” at the Seventh Conference on Polar Meteo- 
rology and Oceanography, Hyannis MA, May 12-16, 2003. 

• Presented a talk “Assimilation of ice thickness information into a sea ice model”. Proceedings of the 
Sixth Conference on Polar Meteorology and Oceanography, Amer. Meteorol. Soc., 15-20 January 
2001, San Diego CA. 

• A monthly averaged gridded data set produced was produced and posted on the Arctic Ocean Model 
Intercomparison (AOMEP) web site and at psc.apl.washington.edu/lindsay/#RGPS. 

• Generated a movie of the model results, a simulation showing the trajectories of 700 cells for a 3- 
month period. It may be found at http://psc.apl.washington.edU/lindsay/#TRAM under animations. 

• We have made substantial progress on formulating the data assimilation techniques we will use to 
assimilate RGPS data and have performed assimilations for the SHEBA year (1997-98). The result- 
ing ice velocities are more accurate (compared to buoys) than other model or data sources (Table 1) 
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4. Future work 


We are currently in the process of preparing a paper describing the data assimilation concepts out- 
lined above and will finish the paper under a pending NASA grant “Sea ice thickness estimates in mod- 
els and observations, the view from space”. We have also submitted a proposal to the Oceans and Ice 
2004 NASA Research Announcement (NRA-04-OES-02) to expand and continue this work. The objec- 
tives of the new proposal include the following four science objectives: 

• Determine the long-term trajectories and evolution of ice from important source regions. 

• Determine the spatial and temporal variability of the mixing and diffusion of the sea ice cover. 

• Enhance our understanding of the spatial structure of the stress and strain through the analysis of 
simulations constrained by assimilation of RGPS data. 

• Enhance our understanding of the errors in wind forcing fields through comparisons of RGPS tra- 
jectories with model output computed with three different wind fields: LABP, NCEP Reanalysis, 
and the European ERA-40 Reanalysis. 

In order to advance these science objectives we will need to address two technical objectives: 

• Improve the accuracy of the Lagrangian model in simulating the ice motion and ice thickness 
through model improvements based on validation with RGPS data, and extend assimilation capabil- 
ities of the model to include ice extent. 

• Provide a gridded and interpolated version of the entire RGPS trajectory data set. It will be the most 
accurate and detailed regularly gridded data set of ice motion available for analysis and for valida- 
tion or assimilation into sea ice models. 

It is hoped that we will be able to continue our work with this exciting new tool for sea ice analysis. 
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Figure 1. Time series of one component of 
the positions of one Lagrangian cell and 
the forces acting on it In the first panel the 
first-cut (model only) trajectory is shown 
in red. The locations of the RGPS trajec- 
tory measurements are shown in blue. A 
curved blue line, derived from the red 
curve, connects them. The green curve is 
the final trajectory that is forced to follow 
the blue line by the addition of a corrective 
force proportional to the difference 
between the green and blue lines. The bot- 
tom panel shows the forces: cyan for wind 
stress, blue for water stress, orange for 
Coriolis force, and red for internal stress. 
Black is the corrective force. 


Figure 2. Map of the model ice velocity 
(green), wind stress (cyan), and corrective 
force (black) for one day. In this case the 
wind stress and ice drift are approximately 
uniform, but the corrective force is mostly 
small except for a block of cells in the center. 
This suggests a local problem with the inter- 
nal stress and a rigid plate. 
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